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Abstract 



Penitentes and suncups are structures formed as snow melts, typically 
high in the mountains. When the snow is dirty, dirt cones and other 
structures can form instead. Building on previous field observations 
and experiments, this work presents a theory of ablation morphologies, 
and the role of surface dirt in determining the structures formed. The 
glaciological literature indicates that sunlight, heating from air, and dirt 
all play a role in the formation of structure on an ablating snow surface. 
The present work formulates a mathematical model for the formation 
of ablation morphologies as a function of measurable parameters. The 
dependence of ablation morphologies on weather conditions and initial 
dirt thickness are studied, focusing on the initial growth of perturbations 
away from a flat surface. We derive a single-parameter expression for 
the melting rate as a function of dirt thickness, which agrees well with a 
set of measurements by Driedger. An interesting result is the prediction 
of a dirt-induced travelling instability for a range of parameters. 
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Penitentes are structures of snow or ice jTJ, which commonly form during the 
summer on glaciers or snow fields at high altitudes (in the Andes and Himalaya). A 
penitente is a column of snow, wider at the base and narrowing to a point at the tip. 
The name "penitente" is a Spanish word meaning "penitent one," and arose because 
a field of penitentes resembles a procession of monks in white robes. Penitentes range 
from one to six meters high with the spacing between columns comparable to their 
height (Figure H). Smaller structures, known as suncups or ablation hollows, can 
be found in lower mountains like the Rockies and the Alps (Figure |2]). Suncups are 
smaller, two cm to half a meter in size. 
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FIG. 1. Photographs of penitentes, from Post and LaChapelle pi, p. 72. Left, penitentes 
on Cerro Negro, Chile. Right, field of penitentes, north slope of Cerro Marmolejo Norte, 
Chile. Note the ice-axe, approximately 80 cm high. In the picture on the left, the snow in 
the hollows has completely melted, exposing the soil underneath. This is a frequent, though 
not universal, feature of penitentes 

The first written record of penitentes comes from Charles Darwin, who observed 
them during his travels in the mountains of Chile [[|. 

"Bold conical hills of red granite rose on each hand; in the valleys there 
were several broad fields of perpetual snow. These frozen masses, during 
the process of thawing, had in some parts been converted into pinnacles 
or columns, which, as they were high and close together, made it difficult 
for the cargo mules to pass. On one of these columns of ice, a frozen horse 
was sticking as on a pedestal, but with its hind legs straight up in the 
air. The animal, I suppose, must have fallen with its head downward into 
a hole, when the snow was continuous, and afterwards the surrounding 
parts must have been removed by the thaw." 



An extensive literature of observations and field experiments has documented 
these ablation morphologies (see for many references). Ablation in this context 
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means removal of snow by melting or sublimation. This contrasts with other pro- 
cesses like wind, avalanches, and rain. There is a consensus about the causes of 
ablation morphologies, although some contradictory claims do exist in the literature. 
For penitentes, bright sunlight and cold, dry weather are apparently required [||, 
while ablation hollows can be formed in three distinct ways, with solar illumination 
important in some settings. For other locations, uniform heating from the air appears 
to be the key effect. The effect of this "sensible" heat transfer (so called because it 
is easily felt with the senses) to the snow depends on whether the snow is clean or 
dirty. Since many readers are likely to be unfamiliar with the glaciological literature, 
I give a brief review here. 

The observational evidence for sunlight-driven formation of penitentes is abun- 
dant. In early work, Matthes Q argued that a variety of ablation forms, from sun 
cups a few inches in size to penitentes many feet deep, are formed by the sun. As 
he pointed out, the formation of the largest penitentes requires strong and prolonged 
solar radiation — the primary reason why penitentes develop only in regions with dry 
summer climates. Matthes also observed that penitentes tilt toward the elevation of 
the midday sun (an observation confirmed by others |],|5|-§||). Such tilting is strong 
evidence that the sun has an important role in the development of structure, because 
the direction of incident radiation provides the symmetry axis in the problem. In 
later work, Lliboutry || observed that incipient penitentes begin as east-west rows. 
Perhaps most important, if the weather is not dominated by direct sunlight — if the 
weather is cloudy H or very windy [{||| — penitentes are observed to decay. In the 
1930s Troll performed an experiment trying to create penitentes J7|. The exact state- 
ment (reported by Lliboutry ||) is "Troll was able to reproduce penitents in Germany 
by shining an electric bulb on fresh snow during a cold, dry night." This supports the 
sunlight mechanism, although to my knowledge no controlled laboratory experiments 
have investigated light-driven structure formation. 

To understand qualitatively how sunlight can cause structure formation, note 
that when light is reflected off the snow, the base of a depression receives more 
reflected light than the neighboring peaks. This drives an instability of the surface 
and the amplitude of a perturbation grows; quantifying this argument will be a main 
goal of this paper. The effects of reflections are considered imporant by several 
observers @,|8|||. This may not be the only required effect. At the high altitudes 
where penitentes commonly form, the air is so cold and dry that sublimation occurs 
instead of melting ||10|| , consistent with the observations that the snow in penitentes is 
quite dry ||||. Lliboutry || claims that the snow in the hollows between penitentes is 
soft and wet, and that temperature variations of 5-10° C exist between the peaks and 
the troughs. This was interpreted to indicate snow sublimating from the peaks and 
melting near the troughs — an effect which accelerates the growth of structure, since 7 
times more heat is required to sublimate a volume of snow than to melt it. Lliboutry 
believes this effect to be crucial for the development of the largest structures, and 
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claims that penitentes only appear at altitudes high enough that sublimation becomes 
important. But other researches report results in disagreement with this 



A different set of observations and experiments have led to a very different claim: 
that solar illumination destroys ablation morphologies, while windy weather promotes 
their growth. Leighly O] argued that heat from air (delivered by wind) leads to the 
formation of ablation polygons (cf. Figure 3). Others state |L2]-|14j] that structures 
do not grow in the presence of direct sunlight. Ashwell and Hannell claim fl5j that 
when the incident solar power is larger than the incident power from heating by wind 
the hollows decay. Detailed observations, along with wind-tunnel experiments, have 
been made by Takahashi and collaborators [0,0 ; they conclude that structures grow 
most rapidly when the air temperature and wind speed are highest | I8| . When the 



weather is warm and cloudy, wind mixes the air so heat is delivered at a steady rate 
to the surface; the higher the temperature and wind speed, the faster the heating. 

Rhodes, Armstrong, and Warren |2| suggested a resolution to this apparent contra- 
diction, which I now summarize. In their view, dirt on the snow surface is the hidden 
variable distinguishing the two cases. Sunlight drives formation of penitentes in clean 
snow because reflection into hollows makes depressions in the snow surface grow. 
Any source of ablation which transfers heat uniformly to the snow surface therefore 
disrupts formation of structure. However, sunlight acts differently on a dirty snow 
surface. Dirt decreases the amount of reflected light, preventing the concentration of 
sunlight in the hollows. This agrees with the Rhodes et. al. observations of suncups 
on Mount Olympus. The researchers noticed that when the snow surface was covered 
by a layer of ash from the eruption of Mount Saint Helens, suncups did not form. 
They scraped away the ash from one patch of snow and observed the formation of 
sun cups on this clean snow surface. 





FIG. 2. Photographs of suncups, from Post and LaChapelle |lj. Left, suncups on the 
Taku Glacier, Coast Mountains, Southeast Alaska, p. 70. Right, deep suncups in Disap- 
pointment Cleaver, Mount Rainier, p. 71. 
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How does dirt affect snow ablation? If the dirt thickness covering the snow is 
sufficiently thick, the dirt forms an insulating layer which slows down the ablation 
rate of the snow pO]. Thus dirt can have different effects, depending on thickness. 
A thin layer of dirt causes faster ablation because reflection is inhibited. However, 
sufficiently thick dirt slows ablation. A large amount of work has looked at how 



debris-covered ice or snow melts [ 10,15,21 1; one typically finds a peak in the ablation 



rate for dirt thickness around 0.5-5 cm. One nice experiment was done by Driedger 



1 22 1, who measured ablation rate as a function of ash thickness on the South Cascade 
Glacier (primarily due to melting). The typical grain sizes of the ash were 0.25 to 
1.0 mm diameter, and the maximum ablation rate occurred for a dirt thickness of 3 
mm. The data from her measurements are shown in Figure || Comparison to these 
data provides a test of my model, as discussed below. 



As pointed out by Ball p3| , small particles of dirt adhere to the snow surface. 



This is true only for sufficiently small dirt particles: the adhesive force on the particle 
must be large compared to the gravitational force |15[]. When adhesion to the snow 



dominates, the pieces of dirt move perpendicular to the snow surface (rather than 
falling straight down) as the snow ablates. Sticky dirt therefore tends to become 
concentrated on the most elevated regions of the surface. (See Figure |].) The con- 
centration of dirt on melting snow can be observed in old snow piles in cities, and is 
illustrated in Figure || This movement of dirt normal to the melting snow surface is 
quantitatively well-documented in the literature (summarized in [Q). For the argu- 
ments here to be correct qualitatively, the dirt need not move completely normal to 
the surface — a component of motion normal to the surface is adequate. 

This mechanism explains dirt-driven structure formation: as the snow ablates, 
dirt becomes concentrated on the more elevated parts of the surface. The thicker dirt 




ablation hollows with dirt collected on the ridges. The structures are reportedly 12 to 18 
inches high, p. 196. Right, dirt cones, approximately 20 to 40 inches high, p. 190. 
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Melting Snow Surface 




FIG. 4. Motion of dirt on a snow surface. A particle adhered to the surface of the snow 
moves normal to the surface (left). When particles follow such "normal trajectories," peaks 
are stable equilibria and valleys unstable equilibria (right). 



forms an insulating layer on the ridges, so they ablate more slowly. The hollows thus 
grow deeper. This concentration of dirt by ablation can, in extreme cases, lead to the 
formation of so-called dirt cones: cones of snow or ice covered by a thick layer of dirt 
T0| , |20|pT| , |24]| . (See Figure |3|.) These structures can become quite large: Swithinbank 



f24| reports a dirt cone in the Himalaya estimated to be 85 m high! Drewry |2T| has 
done detailed experiments on dirt cones. He concludes that the cones ultimately reach 
a steady state, where the motion of dirt toward the center of the cones is balanced 
by the debris sliding down the cone when the slope angle is too large. 

The proposal of Rhodes, Armstrong, and Warren that uniform heating causes 
structure only for dirty snow does not completely resolve the disagreement about 
structure formation. Some observers who advocate uniform-heating driven formation 
of ablation hollows insist that dirt on the snow surface is not required flfllD^IEjITT . 



Indeed, some photographs show ablation hollows in clean snow inside a tunnel or on 
other inverted surfaces, suggesting that neither dirt nor solar illumination are nec- 
essary. How can this be explained? Some have suggested that a regular pattern of 
convection cells leads to the observed polygonal pattern |TT|,[T2 
mate shows this cannot give the correct size structures 
that the structures are formed by turbulent eddies [p"3l ,|T 
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but a simple esti- 
. Another suggestion is 
although Takahashi ]17 



later claimed that the diameters of the hollows are independent of the eddy size. 
Takahashi 0] proposed that the separation of the air boundaray layer as it flows 
over a cusp could produce lower temperatures at the cusp, and therefore lead to 
structure formation. I am not familiar with any further theoretical or experimental 
consideration of the Takahashi proposal; this mechanism for structure formation will 
not be considered in this paper. 

Despite the extensive observations of ablation morphologies, there is a lack of 
mathematical models of their growth | 2~6fl . The goal of this paper is to quantify 



the primary mechanisms discussed above, and characterize the initial stages of the 
instability of a flat surface. We examine how reflection of sunlight and dirt affect 
structure formation. The nonlinear evolution more appropriate to penitentes will be 
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examined in a future paper. 

In this paper we consider sunlight — direct and reflected — the primary source of 
heat which leads to snow ablation. It is well-documented that radiation is the dom- 



inant heat source for ablating snow |27| , |28|| , especially at high altitudes and low lati- 
tudes [0] . The importance of considering reflected light as well as direct illumination 
is supported by the observational evidence. The fraction of light reflected from old 
snow is about 0.5 HlOl^TH ■ Therefore the amount of heat absorbed locally (and cor- 



respondingly the ablation rate) can vary by up to a factor of two for different parts 
of the surface — such a large variation can have important consequences for structure 
formation. Kotlyakov and Lebedeva M made measurements of the albedo on a glacier 
with small penitentes. In a measurement averaged over surface features, ten percent 
more light was absorbed when the sun was high overhead, presumably indicating the 
absorption of reflected light in the structures. 

In the presence of dirt, sensible heating from the air may be important, in addition 
to sunlight. In this paper I focus primarily on the sunlight-dominated case, and 
comment on similarities and differences with sensible heat. Modifying the model to 
include sensible heating is straightforward. 

By forming a quantitative model, we can test whether the effects considered can 
explain the appearance of structure, and describe the morphologies produced. The 
primary goal is to formulate the simplest model which contains the essential physics. 
Ideally the theory would contain no free parameters, that is, all parameters in the 
model can be calculated or measured in experiments. We also discuss which effects 
are left out of the simple model, and estimate how serious are the consequences for 
such omissions. 

The first part of this paper addresses clean snow only. In section [j] we formulate 
a minimal model, and carry out the analysis of the model for small perturbations. 
The linear stability analysis lets us estimate the wavelength of the fastest-growing 
disturbance, and determines the initial size structures that form. This analysis allows 
us to make testable predictions of the conditions necessary to observe structure. 

We then discuss the effect of dirt and re-formulate the model to include dirt in 
section [H[ We compare our model to the field experiment of Driedger and find good 
agreement. Thus reassured that the theory contains the imporant physical effects, 
we show how dirt alters the growth of small perturbations. We show that a thin 
dirt layer suppresses the reflection-driven instability and induces travelling dispersive 
waves on the surface. In the limit of thick dirt, we demonstrate the insulation-driven 
instability expected from the discussion above. 



7 



I. LIGHT REFLECTION ON CLEAN SNOW 



The model for penitente growth we derive here contains simplifying assumptions; 
we hope to capture the essential features while neglecting some effects We will discuss 
the assumptions and their limits of validity as the model is developed. Some of the 
most important simplifications include considering the latent heat to be constant and 
including only first-order, isotropic reflections. We focus on a one-dimensional model 
of penitentes, assuming invariance in the transverse direction, although it is straight- 
forward to generalize these equations to two dimensions or to multiple reflections. 

We consider the height of the snow surface h(x,t), and seek an equation for the 
time evolution of h. 



Heat incident on the surface leads to ablation — the height h decreases as the snow 
melts or sublimates. We assume that ablated snow vanishes into the air or drains, and 
therefore that the flow of water along the surface and re-freezing are not important 
(and similarly that other changes in the nature of the snow are unimportant). This 
model can apply to either melting or sublimation. We use the term "ablation" to 
refer to removal of snow in either way. 

Suppose a point on the surface absorbs a power per unit horizontal area P(x). 
The latent heat required to ablate a unit volume of snow is L. Combining this with 
an effective diffusive smoothing term (see below) gives the evolution equation for the 
surface: 



For clean snow, we assume that L is a constant (independent of x). This is 
true when the surface temperature and humidity are approximately constant. As 
discussed in the introduction, fully developed penitentes may have melting in the 
hollows and sublimation in the tips — a situation which requires L to vary along the 
surface. Indeed, the variation in L might be the essential effect for large structures. 
For small angle structures, that is, amplitude small relative to wavelength, L = 
constant should be a good approximation. Later, we will include spatial variation in 
the effective L due to dirt on the snow surface — see section |I|. 

The second term in equation [I] for the surface height is a simple form of the small- 
scale cutoff: a diffusive term with diffusion constant D. As we will see below, in the 



A. Snow Ablation 
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absence of any smoothing term, the model can produce arbitrarily small structures. 
This is clearly not realistic, because the physics at small scales will cut off the instabil- 
ity. For the qualitative results here, the exact mechanism of the small-scale cutoff is 
not essential; the main point is that there is some minimum size structure which can 
form. A natural small-scale cutoff is the extinction length of sunlight, which defines 
the thickness of the snow layer in which the light scattering takes place [pTfl . Points 
on the snow surface within one extinction length are not optically independent, and 
therefore such nearby points ablate at the same rate. The extinction length depends 
on the density and grain structure of the snow. The typical extinction length [27,|28| 
for old snow (grain radius 1 mm) is of order 1 cm |29[ . We will choose the diffusion 
coefficient so that the characteristic cutoff length is of order the optical extinction 
length. Again, remember that this term in the height equation is a simplified rep- 
resentation of the small-scale physics, and any conclusions which depend sensitively 
on the form of this term should be considered suspect. Note that diffusion of heat 
through the snow might seem another natural form of the small-scale cutoff; however 
the gradients of temperature in the snow are not large enough for thermal diffusion 
to stabilize short wavelengths |30| . 

We note here that recent work by Nodwell and Tiedje |25] considers the scattering 
of light in the snowpack in quantitative detail. They find a range of length scales where 
suncups can form (both a minimum and maximum wavelength), a result which our 
simplified model cannot produce. 



B. Light Reflection 



In this section we describe the reflection of sunlight from the snow surface. We 
assume that the sunlight shines directly down (in the — z direction) and has a uniform 
power per unit length /. The parameter characterizing reflections is the albedo a, 
which denotes the fraction of light reflected. Thus the absorbed power per unit length 
is (1 — a)I. For old snow — called firn — a typical value is a = 0.5 ||f0|| . 

The reflecting properties of snow are different from those of a mirror. Snow looks 
white because it scatters light in many directions, as we would expect for a rough 
surface. Here we treat the light using ray optics, and assume the surface reflects 
isotropically, thus the power is distributed uniformly into n of solid angle outside the 
surface. We approximate that the reflection occurs at the surface of the snow. (As 
mentioned above, the reflection takes place in a layer of order 1 cm thick. We ignore 
this in formulating the reflections, and include its effects schematically through the 
diffusive term.) 

Using these properties, the total amount of light scattered from an interval around 
point X\ to the interval between x and x + dx is 
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71 



d6 dxi, 



(2) 



where dQ is the angle subtended by the surface between x and x + dx (see Figure ||). 

We can find d9 in terms of the shape of the surface. 

dl |p x ds\ Ah — Axh'(x) 
dv = — = — 



V 



V 



Ah 2 + Ax 2 



where we have used h' = dh/dx and 



Ax = Xi — x 

Ah = h(x±) — h(x) 



(3) 
(4) 



and ds is the vector tangent to the surface 

ds = dx(l, h!) 



(5) 



We define the vector p, which points from the point x\ to the point x. From Figure 
[|, we can see that 



p = VAx 2 + Ah 2 



(6) 



Snow Surface 
h(x) 




x' x x + dx 

FIG. 5. Schematic of the ablating snow surface. Scattering from the point x\ to the 
interval between x and x + dx depends on the angle dO. The vector p points from x\ to 
x and the increment dl is normal to p such that d6 = dl/p. The vector n is normal the 
surface at x and ds is the increment along the surface between x and x + dx. 
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To find the total power reflected to point x, we must add up the intensity scattered 
from all points X\\ 

„. . al f dxAAh — h'(x)Ax) , , 

The integrand in this equation is the propagator for light intensity, it describes how 
the intensity is carried from one point to another on the surface; P r {x) is the intensity 
due to a single reflection. To include multiple reflections, we can write the power as 
an integral equation for P. 

, ot f dxiP(xi)(Ah — h'(x)Ax) 
P(x) = ( 1 — a)I 1 



irJ Ax 2 + Ah 2 

This can be written as a power series in a. We will only consider single reflections 
here, which does not introduce a large error when a is small. For old snow, a typical 
value of a ~ 0.5. Including the higher-order correction from multiple reflections may 
be important in determining the precise details of the largest shapes. 

This formula for reflected intensity is not complete, because it neglects the line-of- 
sight constraint. Light cannot scatter from X\ to x if the path of the light ray is blocked 
by another part of the surface. This requirement is a nonlinear constraint which is 
difficult to handle analytically but is straighforward to implement in numerics. We 
typically indicate the constraint schematically, by writing "line of sight" under the 
integral: 



[x) 
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al r dxi(Ah — h'(x)A 
~ 1 Ax 2 + Ah 2 



line of 
sight 



We can also write a necessary (but not sufficient) criterion for the line of sight con- 
straint, when applied to local analysis within one "basin." The two points x and Xi 
are within line of sight of each other when the dot product of the vector normal to the 
surface and the vector p is less than 0: — n • p = Ah — Axh'(x) > 0. (See Figure |5|.) 
Note, however, that this simple criterion will miss intermediate bumps in the surface. 
In other words, the constraint may be satisfied but no reflection occurs between x\ 
and x because the line of sight is blocked by an intervening peak. 



C. Model 



The equations combining reflection and ablation are 
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where we have defined the integral 



1 /■ d Xl (Ah - h'(x)Ax) 
{> nJ»r« Ax 2 + A/i 2 ' 1 ' 

sight 

The intensity of the sun determines a characteristic ablation rate IL" 1 , where L is 
the latent heat per unit volume. Combining this velocity with the diffusion coefficient 
D gives a length 



and time 

DL 2 

t = —- (12) 

The solar constant gives the intensity of solar radiation at the top of the atmosphere 
f28fl I = 1.4 x 10 6 erg cm -2 sec" 1 ; we therefore choose / = 10 6 erg cm -2 sec" 1 as 
the typical value of I under bright sunny conditions. The latent heat depends on 
density. Freshly fallen snow has a density of between 0.05 and 0.2 g cm -3 , while 
older snow that has survived one melt season has a density range of 0.4 to 0.8 g cm -3 
IQySfl. Here we pick an intermediate density of 0.3 g cm" 3 for our estimates. This 



gives a latent heat per unit volume for melting L = 10 9 erg cm 3 and a melting rate 
I/L = 10" 3 cm sec" 1 [31]. We pick D = 2.5 x 10" 5 cm 2 sec -1 , where this choice 



is made so that the most unstable wavelength is 2 cm (see below). In this case the 
length scale I = 0.25 mm. and the time scale t = 25 seconds. 

For sublimating snow, the latent heat is seven times larger, this gives the slower 
melting rate I/L = 1.4 x 10" 3 cm sec -1 , larger length scale i = 1.75 mm and time 
scale t = 1225 seconds. 

We will now perform a perturbation analysis of Equations [| to see how the size 
structures formed compares to the scale £. We have set up the problem so that 
structures will initially form on a scale roughly comparable to £, and expect the 
perturbation analysis to give this result. 



D. Quasi-Linear Regime 



Here we show how an approximate linear analysis of the equations can be per- 
formed. This allows us to derive the dispersion relation, which characterizes when 
the system is stable or unstable. There is a fastest growing mode determined by the 
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competition between reflection and diffusion. The length scale of this mode is related 
to the basic scale i from dimensional analysis above; we determine the pref actor here. 
The results are significant because they describe how the physical parameters affect 
the instability. We will argue that reflection favors structures on scales as small as 
possible. On the other hand, the small-scale cutoff limits the smallest structures 
possible. Therefore we expect the fastest-growing mode to be of order the cutoff size. 

The reflection integral is scale invariant: upon rescaling x and h by the same 
amount the integral 1{x) is unchanged. Thus in the absence of diffusion, there is 
no characteristic scale in the problem. Therefore a shape with aspect ratio 1 — a 
shape with variations in h comparable to variations in x — should have a growth rate 
of order 1 (in the absence of boundary effects). The integral contributes a shape 
factor independent of the amplitude of the shape 5. Therefore the rate of change of 
amplitude 5 is constant. 

To examine shapes with aspect ratio far from one, we start with an aspect-ratio 1 
shape, then transform x — > \x and h — > 5h . When 5 <C A, we find that the integral 
scales with the basic angle 8/X: I — > 5 /XL. Thus for small perturbations, we expect 
a growth rate proportional to the amplitude (5 ~ 5). 

For sufficiently small S/X, we treat the contribution from the reflection integral as 
a numerical factor of order 1. Note that a sinusoidal perturbation is not an eigenshape 
for small amplitudes; we do not know what the actual eigenshapes are. The dominant 
contribution is the scaling with 5, and we neglect the other (slower) dependence 
on position, amplitude, etc. Thus the quasilinear equation for a small-amplitude 
variation in the surface h = 5 sin qx is approximately 

1(x) ~ — 8smqx e w * (13) 
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which gives a dispersion relation 

al 



u = VL q ~ Dq2 (14) 



This argument selects a fastest-growing mode with wavenumber 

al 



2nLD 



(15) 



* = £i?D=* D (16) 



These equations are the dimensional analysis result, with an estimate of the prefactor 
from the scaling argument. Plugging in values of typical parameters given above, we 
find the most unstable wavelength for melting A* = 2n/q* of 2 cm, and characteristic 
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time 4000 sec. In the case of sublimation the wavelength is 14 cm and the time 2 x 10 5 
sec. The choice of the diffusion coefficient is now clear: we chose D to give a most 
unstable wavelength of 2 cm. We have put in diffusion as a simplified representation 
of the small-scale physics, and chosen its value so that the numbers make sense. It 
is important to remember that because of this choice of D, the numbers calculated 
here cannot be considered a prediction of the initial size structures that form. The 
calculation of real interest is how this instability is changed by dirt, as discussed in 
the following section. 

Although it agrees well with simulations of initial growh of perturbations which 
compute the reflected intensity at each point , we must remember that this analysis is 
only quasi-linear because we do not know the eigenfunctions of the reflection integral, 
and superposition does not hold: because the integral is nonlocal, a surface variation 
with two modes of different wavelength cannot be described by the addition of two 
modes with different q. 



II. EFFECTS OF DIRT 



A layer of dirt on the surface of the snow changes its properties. We model both 
the optical and insulating effects of dirt, and fit the theory to melting data measured 
by Driedger p2| . These data allow measurement of a crucial parameter in the model, 



and the good agreement between theory and experiment show that we have captured 
the important effects of dirt. The essential features are that thin dirt speeds ablation, 
because it increases absorption, while thick dirt insulates the snow, slowing ablation. 
This basic behavior leads to the two different regimes of instability ||. 

Dirt looks black because it absorbs light. The presence of dirt effectively decreases 
the surface albedo and therefore increases the fraction of absorbed light. We assume 
light has a probability of being absorbed that is constant per unit thickness of dirt. 
The fraction of light not absorbed by the dirt is e~ s ^ Se |32| , where s is the dirt thickness 



and s e the extinction length in the dirt — typically of order the characteristic dirt 
particle size. Therefore dirt modifies the albedo according to 

a d = ae~ s/s \ (17) 



Note that absorption by the dirt layer is not isotropic — more light will be absorbed 
near grazing incidence, decreasing the reflection even more. The qualitative effect 
of dirt remains the same however, and thus we neglect this anisotropy. Increased 
absorption through a lower effective albedo hastens snow ablation. 

But the dirt also slows ablation. In the presence of an insulating dirt layer, the 
temperature at the surface of the snow is decreased below the ambient temperature, 
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and more heat is required to ablate a given amount of snow. Suppose an amount of 
heat L is necessary to ablate a unit area of clean snow. How much additional heat is 
required in the presence of a dirt layer? At steady state the temperature satisfies 

V 2 T = 0. (18) 



When the radius of curvature of the surface is large compared to the dirt thickness 
(the important limit for growth of perturbations) we can treat the snow surface as 
planar, leading to variations in T in the z direction only. The boundary conditions 
are: At the dirt-air interface (z — 0), the temperature must be equal to the ambient 
temperature. The temperature gradient at the surface due to heat flux into the dirt 
from the air is T'(z = 0) = Pj k, where P is the incident power flux and k the thermal 
conductivity of the dirt. Thus we find that the temperature at the snow surface is 
less than T(z = 0) by an amount AT = Ps/k. An extra amount of heat AQ = CAT 
is needed to raise the snow temperature up to its value in the absence of dirt, where 
C is the heat capacity of the snow. Thus the effective latent heat for a dirt thickness 
s is 

L d = L + ^ (19) 



Both L and C depend on the ambient temperature T. However, the dependence is 
sufficiently weak that we can neglect it. 

Combining these two effects we find that the snow ablation velocity for a flat 
surface covered with dirt is 

m(s) = ^g(s) (20) 

where g is a dimensionless function of the dirt thickness. In this model, 

1 - ae~ s/Se . , 

g(s) = 7 (21) 

yK ' l + 7s(l -ae~ s /^) y ' 



where we have defined the dimensionless measure of the insulating value of dirt: 

i = "ST < 22 > 

The non-monotonic behavior of this curve — positive slope for small s and negative 
for large s — is the important qualitative result. Note that in the absence of dirt the 
ablation rate is as expected: 

m( s = 0) = -(1 - a) (23) 

Lj 
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— Theory 

□ Driedger 1981 
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FIG. 6. A plot of the relative ablation rate m/m{s = 0) versus dirt thickness. The 
points are the data measured by Driedger |^2[ . The solid curve is a one-parameter fit to 
Equation |2l], yielding the fitted 7 = 0.047. Note the fastest ablation occurs for dimensionless 
s ~ 3. We picked s e = 1 mm from Driedger's measurement of the dirt particle size, and 
albedo a = 0.5 from other measurements [27]. We can also estimate the parameter 7 (see 
text). The estimate gives 7 within a factor of 2 of the value obtained from this fit. 



A fit to the data of Driedger |22| is shown in Figure g. Driedger measured melting 
rates of a flat surface for different dirt thickness. The plot shows m/m(s = 0) versus 
dimensionless dirt thickness. We picked s e = 1 mm from Driedger's measurement 
of the dirt particle size. Fitting the data to Equation [H] allows us to determine the 
dimensionless insulation coefficient 7. For a ~ 0.5 (from other measurements p7| ) 
the fit gives 7 = 0.047. We can also estimate the parameter 7 using other data (see 
below). The estimate gives 7 close to the value obtained from this fit. 

This model and the experiment of Driedger are in the regime where solar radiation 
is the dominant heat source. The discussion of Rhodes et. al. || points out that 
the ablation curve changes when sensible heating is important. In fact, if radiation 
is negligible the curve will monotonically decrease as the dirt thickness increases, 
because light absorption effects disappear in this limit. It is straightforward to adjust 
the model to include other sources of heating. Measurements of the type Driedger 
performed, compared to the type of model presented here, could in principle give 
information on the relative importance of radiant and sensible heating. 
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A. Dynamics of Dirt 



As the snow surface ablates, the dirt layer on it moves (Figure [|). We assume 
the particles are sufficiently small that the snow moves purely normal to the surface. 
The sideways (x-direction) velocity of a piece of dirt is 

v = -hti (24) 



The thickness of the dirt s(x) must obey a conservation equation s+V- (vs) = 0, since 
we assume dirt is neither deposited on nor removed from the surface. The evolution 
equation for the thickness of dirt is thus 

| = -£<«.) = OMsY (25) 



When the surface of the snow is flat {h! = 0) the velocity of the dirt v = 0. Thus the 
tops of peaks and the bottoms of valleys are equilibrium points. The peaks are stable 
equilibria where dirt becomes concentrated, while valleys are unstable (Figure 



B. Model 



We now rewrite the model equations incorporating dirt. We have equations for 
the height of the surface h, the dirt thickness s, and the incident power P. 



^ = Pp . 



(hh's)' 



(27) 



The only sources of heat flux P we will consider are direct and reflected radiation. 
P(x) _ a/ I ae- s ^I f dxi (Ah-h'(x)Ax) 



L L ixL J line ° f Ax 2 + Ah 2 

sight 



We use the same reference ablation rate as in Section |: I/L = 10~ 3 cm sec -1 . 
However, the presence of dirt introduces a new length scale in the problem: the length 
scale for light absorption by the dirt. We choose to nondimensionalize in terms of this 
length, since the physically important regimes of thin and thick dirt are measured 
relative to this thickness. When Driedger measured diameters of ash particles on 
a glacier, 90 percent of the particles had diameters between 0.25 and 1.0 mm |[22|| . 
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We therefore choose s e — I mm as the order of magnitude extinction length for dirt 
absorption; this choice is supported by the good fit to the data. 

The dimensionless timescale comes from combining the ablation rate and length 
scale: id = Ls e /I = 100 seconds. This is the time for a depth s e of snow to melt 
in bright sun. Glacial debris has kk2x 10 4 erg cm -1 sec -1 °K This allows 



us to estimate the dimensionless parameter 7 = s e CI /(Lk) = 0.03. Note that the 
thermal conductivity and the specific heat depend on the density, wetness, etc. The 
fit to Driedger's data (Figure gives a value of 7 ~ 0.047, somewhat larger than 
this estimate. We interpret this as a measurement of the dirt thermal conductivity 
k, and therefore use the implied value k = 1.3 x 10 4 erg cm _1 sec _1 °K. The nondi- 
mensionalized diffusion constant is Did/ s 2 e = 0.25. 

For sublimation the time scale id ~ 700 seconds and the dimensionless diffusion 
constant Did/ s 2 e ~ 1.75; the dimensionless parameter 7 similarly decreases by a factor 
of 7. 

The nondimensionalized equations are 

h = - - P - +DV 2 h (29) 
1 + -yPs v ' 

s=(hh's)' (30) 
P-r(l- ae- s ) + ^-^ f dXl (A/ * ~ h '^ Ax ^ 

sight 

The dimensionless control parameters are r, the solar light intensity; and s, the initial 
dirt thickness. Here we have introduced the parameter r: 

r = L 10 3 sec/cm (32) 

Lj 

to examine the effects of varying the light intensity away from the typical value. 



C. Linear Analysis 



Here we analyze the stability of equations (P9(— |3T|), including effects of dirt. There 
are two important regimes: when the initial dirt thickness is small compared to s e , 
the dirt acts to modify the reflection-driven instability. We find that the instability is 
suppressed by the absorption of the dirt layer, exponentially in the dirt thickness. In 
this regime, dirt can also a induce travelling, dispersive instability of the snow surface. 
Qualitatively, this dispersion arises from the coupling of dirt motion to absorption. 
Dirt migrates to the highest point on the surface — but then the thicker dirt increases 
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FIG. 7. The symmetric (left) and antisymmetric (right) modes of dirt modulation. 
The antisymmetric mode is the physically important one because the symmetric mode is 
unstable. 

the ablation of that peak, and it ablates until it is no longer a local maximum. The 
existence of these waves is an experimentally testable prediction which has not, to 
my knowledge, been discussed before. 

The other limit is when the dirt thickness is large compared to s e . The effective 
albedo ae~ s — > 0. Therefore the dirt instability is independent of light reflections; 
the "light" therefore acts simply as a source of heat. The instability is driven by 
dirt insulating the snow. The characteristic length and time scale of the instability 
depends only on the thermal properties of the dirt. Within this insulation-dominated 
regime, the behavior of the instability depends on whether s « 1/ (jr) or s > 
1 / {p(r) — see below. Thus there are three different regimes of behavior, depending the 
dirt thickness. 

As mentioned above, under different weather conditions uniform heating from the 
air may be more important than radiant heating. In this case any amount of dirt slows 
ablation of the snow |2|], and the insulation-driven instability is the only one possible. 
This can be included in the model by removing the dirt- dependent absorption of light. 

We will perform a linear perturbation analysis: we assume that variations of the 
dirt thickness As are always small. However, the initial uniform dirt thickness s D may 
be large or small relative to s e ; this initial thickness determines the limit of instability. 



D. Thin Dirt Limit 

Here we consider the limit s Q <C s e , meaning the initial uniform dirt thickness is 
small compared to the extinction length. 

There are in general two modes of dirt modulation (Figure |7|): the symmetric 
mode with constant thickness and the antisymmetric mode with As = 2e cos qx. The 
symmetric mode, because it has constant thickness, it simpler to analyze. Note that 
constant dirt thickness is unstable: any modulation in the dirt thickness tends to 
grow. 
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1. Symmetric Mode 



Because the symmetric mode has constant thickness, it insulates the snow surface 
uniformly. Therefore, no thick dirt instability can arise from the symmetric mode. 
But the symmetric mode affects the reflection-driven instability. We look for solutions 
of the form 

h = -rat + 5e ut cos qx (33) 



where m(s) is the ablation rate of a flat surface covered with dirt, calculated above. 
The dirt thickness s Q = constant. We expand the equations to first order in S. The 
resulting dispersion relation is 

^ = are~ s °q ^ 2 

tt(1 + (1 - a)^rs ) 2 

Compare this to the clean snow dispersion relation, Equation [14|. The first term 
(proportional to q) contains the factor e~ So . This term decreases exponentially with 
increasing dirt thickness. For s Q much larger than one, this term is so small that the 
instability practically does not exist. The factor (1 + (1 — a)7rs D ) 2 in the dispersion 
relation results from uniform insulation by the dirt layer. 

The most unstable mode q„ is 

®re- s ° 

q * 2vr J D(l + (l-«) 7 r So ) 2 1 ' 

oo, = qlD (36) 



Figure [8| shows how dirt cuts off the instability, with fixed light intensity r = 1. When 
s Q <C 1, the wavelength is close to the wavelength in the absence of dirt. However, 
the absorption of light by dirt becomes important for s > 0.1 and the wavelength 
increases exponentially. As the wavelength increases, the growth rate of the instability 
decreases, and the instability becomes less readily observed. 



2. Antisymmetric Mode 



The antisymmetric mode involves variations in the thickness of the dirt. We must 
solve for the coupling between snow ablation and dirt motion. The solution is of the 
form 

h = -mt + 5e wt cos qx (37) 
s = s + 2ee ut coaqx (38) 
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FIG. 8. Wavelength of the Reflection-Driven Instability. The wavelength is normalized 
to the most unstable wavelength of clean snow. When s„<l, the wavelength is close to 
the wavelength in the absence of dirt. However, the absorption of light by dirt becomes 
important for s Q > 0.1 and the wavelength increases rapidly. The plot is for fixed solar 
intensity r = 1, a typical value. Parameter values for this plot are as discussed in the text: 
a = 0.5, D = 0.25, 7 = 0.047. 



where s is the uniform dirt thickness at t — 0. Upon linearization, Equation |30] for 
the motion of dirt relates the perturbation amplitudes 



e ms q 
6 = 



2uj 



(39) 



The dispersion relation, to second order in q, is 



UJ 



[1± ff\ 



are °q 



2tt(1 + 7r's D ) 



(40) 



where / is, defining w = 1/(1 + (1 — a) , ~frs ) and recalling m = (1 — ae s °)rw is the 
dimensionless melting rate as a function of dirt thickness, 



/ = (are-^ir-WYil + 



4s m(m 2 7 — are So w 2 ) 
(are~ So Tr~ 1 w 2 ) 2 



(41) 



In the limit s a — ► 0, this dispersion relation is identical to the symmetric mode. 
However, for increasing dirt thickness it contains effects from the dirt modulation. 
The term / can be negative, leading to an oscillatory component to u. Thus dirt can 
cause the instability to travel on the snow surface, in a region of phase space shown 
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FIG. 9. Travelling Instability in the Reflection-Driven Instability. Under the line, there 
is an imaginary part of u>, showing the regime where travelling waves exist. The dirt 
thickness is normalized so s Q = 1 corresponds to one extinction length; similarly, r = 1 is 
a typical intensity of sunlight. For the typical solar brightness r = 1, any dirt thickness 
s > 0.008 will show a travelling instability; therefore most dirty snow surfaces should show 
this behavior. Parameter values for this plot are as discussed in the text: a = 0.5, D = 0.25, 
7 = 0.047. 



in Figure f|. For the typical solar brightness r = 1, any dirt thickness s Q > 0.008 
will induce travelling; therefore, most dirty snow surfaces should show this behavior. 
Qualitatively, this arises from the coupling of dirt motion to absorption. Dirt migrates 
to the highest point on the surface — but then the thicker dirt increases the ablation 
of that peak, and it ablates until it is no longer a local maximum. The positive and 
negative roots in the dispersion relation correspond to left and right moving modes. 
The existence of these travelling instabilities is an experimentally testable prediction. 

Note that the equation is not well-behaved for / = 0. When / = the terms 
in the equation coupling motion of dirt to ablation vanish; the dispersion relation 
reduces to the expression for the symmetric mode above. 

When / is negative, we can find the fastest growing wavelength by looking at the 
real part of u: 

ae~ So r 

q * = 27r J D(l + (l-«) 7 rs ) 2 ( 2 > 

^ = jqt (43) 
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E. Thick Dirt Limit 



The equations are considerably simplified in the limit of thick dirt s ^> 1. The 
effective albedo ae~ So — > 0. Therefore the dirt instability is independent of any 
reflections; the quasi-linearized equations are truly linear in this limit. The thick-dirt 
instability is driven purely by dirt motion coupled to slower ablation under a thicker 
dirt layer. This instability is the linear precursor to the dirt cones of Figure |3|. 

Note that if light is not an important source of heat, the "thick dirt limit" is 
actually valid for all dirt thicknesses. 
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FIG. 10. Most unstable wavelength A* versus dirt thickness s , with typical heat flux 
r = 1. The wavelength is normalized to the most unstable wavelength of clean snow. 
Comparing this figure to the thin-dirt instability, we see that when s Q > 1 the wavelength 
will initially decrease, then increase beyond s Q = 20. The growth rate of this instability will 
be greatest where the wavelength is smallest. 



Replacing ae So — > 0, the symmetric mode disappears. The background ablation 
rate m = r/(l + 7rs Q ). The dispersion relation is, to second order in q, 



y/-ys m 3 



D 2 

-q q 

' 2 H 



(44) 



Here no imaginary component to the dispersion relation is present; it is a straight- 
forward linear instability with one growing mode. The most unstable wavenumber 
is 



7r 3 s 



2D\ (l + 7 rs D ) ; 



(45) 
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For a fixed value of the heat input r, the most unstable wavelength scales differently 
at small and large s Q : 

A* ~ s- 1/2 for s < l/( 7 r) (46) 
~ s Q for s > l/( 7 r) (47) 

The location of the minimum wavelength is determined by the dimensionless param- 
eter 7, which represents how well the snow insulates per unit thickness. Therefore, 
even for optically thick dirt s Q ^> 1 there is a change in the behavior, depending on 
the value of s a compared to the insulation parameter. Since typically = 0.05, these 
limits are consistent. 

There is an optimal s D ~ 1/(7?") ~ 2 cm where the wavelength is smallest. Figure 
r0| illustrates this: it shows the unstable wavelength vs. dirt thickness for the typical 
r = 1, with the optimal s a ~ 20 ~ 2 cm. Comparing this figure to the thin-dirt 
instability, we see that when s Q > 1 the wavelength will initially decrease, then 
increase beyond s Q = 20. The growth rate of this instability will be greatest where 
the wavelength is smallest. 



III. DISCUSSION: COMPARISON TO EXPERIMENT 



This paper has presented work on a simple theory to describe the initial formation 
ablation structures such as suncups, penitentes and dirt cones. We have tried to 
make the model as simple as possible while including the essential physics. As we 
have shown, most parameters in the equations can be calculated or measured in 
experiments, allowing predictions with no free parameters. The exception is the 
effective diffusion coefficient D, which we estimate using the value for light diffusion. 
However, we have not realistically treated the small-scale scattering of light in these 
schematic results. 

At this point, the only quantitative comparison between this model and experi- 
ment is the prediction of ablation rate of a flat snow surface, compared with the data 
of Driedger in Figure |6|. This measurement allows us to extract the dimensionless 
constant governing dirt insulation. The good agreement indicates we have captured 
the important effects of dirt. 

The linear stability analysis of the equations shows the two types of instability 
described in the literature. The model predicts the dependence of the most unstable 
wavelength and characteristic growth rate on the experimental control parameters, 
predictions which could be tested. We argue that for little or no surface dirt, light 
reflection drives the instability. This instability is exponentially suppressed by a dirt 
layer, consistent with field observations. We predict travelling modes induced by a 
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modulated dirt layer in this regime. The existence of such travelling modes is an 
experimentally testable new phenomenon. 

In the presence of a thick layer of dirt, our analysis finds the insulation- driven in- 
stability, as expected. Here we predict an optimal dirt thickness where the instability 
is most easily observed, which depends on the thermal properties of the dirt. 

The visually striking structures in the field are the larger structures: penitentes 
and dirt cones. Understanding the nonlinear regime of the model presented here is 
therefore of interest, and will be the subject of a future paper. The scale of both pen- 
itentes and dirt cones is typically larger than the size of smaller-amplitude structures. 
One way to explain this, which has been suggested from observations is that 

large structures grow at the expense of small ones. Such coarsening behavior is also 
apparent in preliminary work on the nonlinear regime of the model presented here. 

The most obvious problem with the results here is that we have considered vari- 
ation of the surface height in only one direction. Checking whether the results are 
the same for a realistic 2D surface is a necessary extension of this work. A better 
understanding of the small-scale cutoff is also important. In particular, we need to 
understand how using different representations of the short-scale physics affect the 
numerical predictions (of the fastest-growing wavelength, for example). 

Because the model here is simplified, we have left out some physical effects which 
may be important in the experiment. Our treatment of light reflection considered 
single reflections only, which may be a bad approximation with the albedo is close to 
1 (large amount reflected). In the field, the sun of course is not always high overhead — 
the variation of the angle of incident light over the course of the day might change the 
shapes. Other possibly important effects which can occur in field situations include 
other sources of heat transfer to the surface, gravity, and the deposition/removal of 
dirt. Better comparison with lab or field experiments should indicate which of these 
effects are most important to include. 
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